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Preface 

Within  the  past  5-7  years,  the  ability  to  model  the 
supersonic/hypersonic  viscous  flows  using  Computational  Fluid 
Dynamics  (CFD>  techniques  has  advanced  considerably.  An 
algorithm  due  to  Dr.  James  Thomas  of  NASA  Langley  Research 
Center,  Hampton,  Virginia,  was  used  in  the  current  study  for 
performing  accurate  and  efficient  flow  field  calculations 
over  blunt  bodies.  For  high  speed  flow,  the  physics  of  the 
flow  can  be  modeled  using  the  Approximate  Navier-Stokes  (ANS) 
equations.  A  relaxation  type  algorithm  is  developed  in  the 
present  work  for  solving  the  ANS  equations  and  solutions  are 
obtained  for  some  model  problems. 

I  would  like  to  thank  Dr.  A.  Halim  for  supervising  this 
study  and  for  his  technical  assistance  and  aid  with  this 
research.  I  wish  to  thank  Mr  C.  Rumsey  of  the  NASA  Langley 
Research  Center  for  checking  the  Thomas  code  on  the  NASA  Ames 
Cray  II.  I  also  wish  to  thank  Capt  Carlson,  for  writing  the 
block  tridiagonal  matrix  solver  used  in  the  ANS  algorithm. 

During  this  study  I  used  the  Cyber  and  Cray  XMP 
computers  at  WPAFB,  OH  for  the  computational  results 
presented  in  this  thesis.  I  also  used  an  Amiga  personal 
computer,  with  Textcraft  word  processing  software  to  write 
this  thesis.  The  thesis  was  printed  on  a  Juki  5510  printer. 
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Thermal  deformations  induced  by  aerodynamic  heating 
high  speed  vehicles  are  an  important  concern  in  design. 
Aerodynamic  heating  may  have  a  significant  effect  on 


performance  of  the  vehicle,  and  effective  techniques  for 
predicting  the  heat  transfer  and  flow  properties  are 


portions  of  the  flow  domain  as  in  the  case  when  the  viscous 
layer  and  the  shock  layer  are  completely  merged,  the 
classical  ways  of  using  the  boundary-layer  interaction 


flight  vehicles  in  high  altitude  regimes  where  nonequilibrium 
chemistry  effects  are  important.  The  perfect  gas  assumption 
has  been  used  for  the  present  investigation.  Different 
approaches  are  available  for  solving  such  problems.  The  NS 
formulations  lead  to  a  nonlinear  system  of  equations.  Using 
an  explicit  algorithm  to  solve  such  problems  results  in  the 
requirement  of  very  small  time-steps  in  order  to  satisfy  the 
stablilty  bounds.  Therefore,  many  iterations  and  large 
computer  times  are  required  to  reach  the  steady  state.  To 
remove  the  time-step  restriction,  fully  implicit  methods  have 
been  investigated.  The  fully  implicit  code  used  in  this 
study  for  solving  high  speed  viscous  flows  is  described  in 
reference  (19).  The  algorithm  uses  implicit  upwind  finite 
volume  flux  vector  splitting  on  the  inviscid  terms  and  second 
order  central  differencing  on  the  viscous  terms.  An 
attractive  feature  of  upwind  flux  vector  splitting  schemes  is 
that  they  are  naturally  dissipative,  and  artifical  viscosity 
terms  are  no  longer  required  to  overcome  instabilities  in 
regions  of  strong  gradients.  Use  of  the  conservation  law 
form  of  the  governing  equations  in  this  algorithm  also  allows 
shocks  to  be  captured,  and  eliminates  the  need  to  use 
shock-fitting  techniques  to  obtain  the  location  and  strength 
of  the  shocks  in  the  flow  (7:1). 

The  inviscid  flux  vectors  are  split  into  forward  and 
backward  flux  vectors  by  splitting  the  eigenvalues  of  the 
Jacobian  matrix  into  positive  and  negative  diagonal 


eigenvalue  matrices.  These  flux  vectors  are  then  differenced 
using  the  appropriate  upwind  or  downwind  scheme  (7:1-2). 

The  algorithm  will  be  used  to  compute  two  hypersonic 
viscous  flows  over  blunt  bodies  using  the  full  NS  equations. 
Results  are  obtained  first  for  the  blunt  body  flow 
surrounding  a  circular  cylinder  which  was  studied  by 
Tannehill  et .  al .  (22).  The  freestream  conditions  chosen  for 
these  computations  are  f'^^4.6,  Re^lOOOO,  Pr*0.72, 
Pej*14.93N/m=,  T-167K  and  X*l*4,  with  a  cylinder  diameter  (D) 
of  0.3048  m  and  a  wall  temperature  of  556K.  Comparison  with 
available  experimental  and  numerical  data  will  be  given. 
Solutions  also  are  obtained  for  supersonic  viscous  flows  past 
a  circular  wedge  geometry.  Qualitative  comparisons  with 
solutions  of  Bey  et .  al.  (10)  using  other  algorithms  will  be 
given.  Bey  et.  al.  used  a  finite  element  algorithm  to 
predict  the  inviscid  flow  field  around  the  same  geometry. 

For  high  speed  flows,  the  Reynolds  number  is  usually 
very  large  and  consequently  the  thickness  of  the  viscous 
layer  is  very  small.  For  problems  of  this  class,  the 
solution  of  the  full  NS  equations  is  unnecessary  and  time 
consuming.  Different  levels  of  approximations  for  the  NS 
equations  became  the  debate  within  the  CFD  community. 

Boundary-Layer  equations,  are  a  simplified  form  of  the 
NS  equations  where  the  pressure  gradient  normal  to  the  body 
surface  is  neglected.  Viscous  terms  with  derivatives  in  the 
direction  tangent  to  the  body  are  also  neglected  after  an 
order  of  magnitude  analysis  is  performed  on  the  NS  equations. 


The  order  of  magnitude  analysis  is  accomplished  by  assuming 
the  velocity  component  normal  to  the  body  is  small  compared 
to  the  velocity  component  tangent  to  the  body.  Although  the 
BL  equations  require  much  less  computational  effort  than  the 
NS  equations,  the  BL  equations  are  limited  to  the  type  of 
flows  they  can  physically  model  accurately.  Thus  an 
intermediate  set  of  governing  equations  were  developed. 

Approximate  Navier-Stokes  equations(  ANS ),  are  also 
simplified  NS  equations  and  like  the  BL  equations  neglect 
viscous  terms  with  derivatives  tangent  to  the  body  but  retain 
all  other  terms.  An  advantage  of  retaining  terms  which  were 
neglected  in  the  BL  equations,  is  that  separated  and  reverse 
flow  regions  can  now  be  computed  (8:421).  The  ANS  equations 
contain  all  the  Euler  equation  terms,  so  that  the  interaction 
between  the  viscous  and  inviscid  regions  of  the  flow  are 
automatically  tak  *n  into  account  (16:1).  The  ANS  equations 
like  the  BL  equations  require  less  computational  effort  than 
the  NS  equations,  and  derivatives  in  the  direction  tangent  to 
the  body  can  be  approximated  using  a  technique  that  marches 
the  solution  in  the  direction  tangent  to  the  body. 

Upwind  relaxation  algorithms  became  an  attractive  tool 
to  solve  the  ANS  equations  and  are  used  by  many  researchers 
(25).  In  the  present  effort,  an  upwind  pseudo-time 
relaxation  algorithm  that  globally  sweeps  over  the  flow  field 
in  the  direction  tangent  to  the  body,  is  used  to  calculate 
steady  state  viscous  flow  solutions.  This  algorithm  is 


similar  to  the  work  of  Thomas  and  Walter  (25),  using  a  first 
order  upwind  scheme  in  the  direction  tangent  to  the  body  and 
a  second  order  central  scheme  on  all  terms  in  the  direction 
normal  to  the  body. 

The  first  model  problem  solved  using  the  developed 
algorithm  is  the  Couette  flow  problem.  This  problem  was 
chosen  since  it  has  an  exact  solution.  Comparison  of  the 
present  results  with  the  exact  solution  provides  a  test  of 
the  validity  and  correctness  of  the  present  analysis  and 
solution  procedure.  Solutions  are  then  obtained  for 
supersonic  flow  over  a  flat  plate  and  comparison  with  other 
solutions  will  be  given. 
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II  Analysis 


Governing  Equations 

The  time-dependent  compressible  Navier-Stokes  (NS) 
equations  express  the  conservation  of  mass,  momentum,  and 
energy  for  an  ideal  gas  in  the  absence  of  external  forces 
(25:1).  The  momentum  equations  use  Newton's  viscosity  law, 
and  the  energy  equation  uses  Fourier's  conduction  law  so  that 
the  NS  equations  can  be  written  in  the  non-dimensional  form 
in  conservation  law  form  and  Cartesian  coordinates  as  follows 

continuity: 


p  u  )*.+(  p  v  )„*0 


x  momentum: 


<  P  u  >*  +  <  p  u*+p-  >„+<  p  uv-  T )x=0 


y  momentum: 


<  pv  )*+(  p uv-  r xy  )K+(  /5va+p-ryy)>,=o 


energy: 


<  e«  )*  +  [  (  e*+p  )u-u  r»„-v  r Ky+q«  )«♦ 


l  <  e%+p  )v-u  r *y-v  ryy+qy  1^  =  0 


where 


q*  -  -fJLCu/l  RePr(  7-1 )  )  (5) 

q„  *  -flCy/l  RePr(  / -1 )  1  (6) 

rxx  -  <  2ux-vy  )2/X /3Re  (7) 

=  (  Uj,  +  vx  )fi  /Re  (8) 

=  (  2v„-ux  )2/0./3Re  (9) 

Re  -  p~c~Lm/  fl~  (  10  ) 

Pr  -  c„m (11) 

c*  -  7r/<  7  -1  >  (12) 

/X  -  1 . 458X10"*kg/(  ms(  K  )*'*  )T*3'*/(  T~+110.4K  )  (13) 


the  definitions  of  variables  are  given  in  the  List  of  Symbols 
and  references  (8:189)  and  (  8:480-482  ).  The  perfect  gas 
equation  is  used  for  closure 

p  *  (/-DPe  (  14  ) 


t  =  /  p /p 


(  15  ) 


Equations  (1,2,3  and  4)  can  be  combined  in  vector  form  and 
written  in  the  strong  conservative  form  as  follows 


Q«  +  (  B-Ew  )x  +  ( Pi_Fv  “  0 


where 


«  =  I  p,  p  u,  pv,  e* 


E  =  l  p  u,  p  ua+p,  p  uv ,  (efc+p)u]T 


Ev  ”  [0,  1  MM  ,  7T "  y  p  U^"m«+V7T«  y  — q*  ] 


F  i  =  [  pv,  p  uv,  pva+ p,  <e%+p)v]T 


Fv  —  [  0,  m  y  t  ly  y  ,  U  1”  n  ^  V  ^"yy  1  ^ 


(  16  ) 


(  17  ) 


(  18  ) 


(  19  ) 


(  20  ) 


(  21  ) 


.  Equations  (1-21)  have  been  nondimens lonai ized  using  the 
following  identities 


t-f/L'/cT 

CD 


u  =  u**/c* 

CD 

P  *  P-'Pa; 

P‘P 


x  =xm/L* 


v  =  v*Vc:« 


T-TVr 

CD 


L"*l 


Y  =  y’*/L* 
e  =  e**/c”* 

CD 

u  -  ^ 

y  rt* 

CD  CD 


i  22  ) 


.  The  highest  order  terms  come  from  the  viscous  forces  which 
are  second  order,  but  there  are  first  order  convective  terms 


Tj  =  7)  (  x  ,  y  ) 


(  24  ) 


where  ^  and  7]  are  defined  in  the  List  of  Symbols.  When  the 
above  transformation  is  applied  to  equation  (16)  the 
following  equation  results 


+  F7* 


(  25  ) 


where 


J  ~  ^  7)  yr'TJ  **  1 


Q’  =  Q/J 


B’  =  (  E-Ew  )+  §x<Fi-Fv)  )/J 


F*  =  [  E-E~  )+T7y<  Fi-F~  )  }/J 


(  26  ) 


(  27  ) 


C  28  ) 


(  29  ) 


J  is  the  transformation  Jacobian.  See  Appendix  A,  for 
details  on  the  transformation  of  the  governing  equations  into 
generalized  coordinates. 

With  the  NS  equations  written  in  generalized 
coordinates,  it  is  easy  to  write  an  algorithm  that  can 
compute  solutions  for  a  variety  of  geometries.  The  method 
described  in  reference  (19)  uses  a  finite-volume  method 
describing  the  balance  of  mass,  momentum,  and  energy  over  an 
arbitrary  control  volume.  The  vectors  ^X/J,  ^y /J ,  7^/J  and 


7]y/J  represent  directed  areas  of  cell  Interfaces  in  the  ^ 
and  7]  directions  in  computational  space.  The  Jacobian 
represents  the  inverse  of  the  cell  volume,  and  the  elements 
p  u/J  and  pv/J  are  the  mass  fluxes  crossing  the  cell 
interfaces  in  the  ^  and  T)  directions.  These  fluxes  are 
then  split  into  forward  and  backward  contributions  to  the 
flow  crossing  the  cell  interfaces,  which  allows  the  algorithm 
to  accurately  resolve  seperated  and  subsonic  flow  regions 
(  6:2  ). 


Approximate  Navler-Stokes  Equations 

The  Approximate  Navier-Stokes  ( ANS )  equations  are  formed 
from  the  NS  equations,  by  doing  an  order  of  magnitude 
analysis  on  the  equations  after  assuming  the  velocity 
component  normal  to  the  body  is  much  less  than  the  velocity 
component  tangent  to  the  body.  After  doing  this  order  of 
magnitude  analysis,  viscous  terms  with  derivatives  in  the 
direction  tangent  to  the  body  are  considered  small  compared 
to  the  viscous  terms  with  derivatives  normal  to  the  body. 

Once  this  is  done  equation  < 25 >  can  be  written  as 

Q'c  +  Ef"  +  =  0  (  30  ) 

where 


E 


(  )/J 


(  31  ) 


F*  =  t  77 *<  E-E"v  >+7^<Fi-F"v  )  ]/J 


(  32  ) 


E‘V  =  to,  t*»k(  rH«Xf  u  r''--+v  r ''-x-q'"  ]' 


(  33  ) 


F”v  *  to,  r"«x,  r  "xx,  u  r  "-x+v  r“xx-q'x  )' 


(  34  ) 


r*'--  =  <  27^mU?  -  7)yV^  )2yU, /3Re 


(  35  ) 


Tr"«x  =  (  7)yu7  +  7)i,Vj  )/X/Re 


(  36  ) 


77  "  xx  -  (27]yv?-7),u7  )2/X/3Re 


q'«  -  RePr<  / -1 )  ] 


(  37  ) 


(  38  ) 


q ' x  =  -  ./(  RePr(  T-l  )  ) 


(  39  ) 


.  It  has  been  demonstrated  that  these  ANS  equations  are 
applicable  to  a  wide  class  of  problems  (25:2). 

The  main  feature  of  the  ANS  equations  is  the  fact  that 
the  highest  derivative  in  the  streamwise  direction  is  less  by 
one  than  those  for  NS  equations.  That  means  the  ANS 
equations  are  parabolic  for  the  viscous  flow  region  but  they 
are  also  hyperbolic  for  the  inviscid  flow  if  it  is 
supersonic.  So  the  equations  are  actually  a  mixed  set  of 
hyperbolic-  parablolic  equations.  That  is  the  reason  the 
abbreviation  ANS  was  preferred  over  Parabolized  Navier-Stokes 


( PNS  )  or  Partially  Parabolized  Navier-Stokes  ( PPNS  ) 
abbreviations.  The  ANS  equations  have  recently  become 
famous,  because  they  can  predict  complex  three  dimensional, 
steady,  supersonic  viscous  flow  fields  in  an  efficient 
method.  The  efficiency  for  three  dimensional  solutions  came 
from  the  fact  that  the  equations  can  be  solved  using  a 
space-marching  finite  difference  technique  as  opposed  to  the 
time  marching  technique  normally  employed  for  NS  equations. 
To  be  able  to  use  the  ANS  equations,  the  solution  at  a 
certain  plane  or  location  in  the  streamwise  direction  has  to 
be  known  as  opposed  to  the  NS  equations.  The  details  of  the 
numerics  are  shown  in  Chapter  III. 
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III  Numerical  Solution  of  the  Equations 


The  important  aspects  concerning  the  numerical  solution 
of  the  governing  equations  are  discussed  in  this  chapter. 

The  finite  difference  form  of  the  equations  are  presented, 
together  with  the  int ial izat ions  used  for  all  the  variables. 
The  procedure  for  solving  the  resulting  algbraic  simultaneous 
equations  is  described  and  some  details  about  the 
implementation  of  the  boundary  conditions  are  discussed. 
Finally,  the  pertinent  numerical  details,  such  as  the  choice 
of  suitable  time  steps,  the  definition  of  convergence,  etc., 
are  presented. 

Fully  Implicit  Scheme 

After  flux  splitting  the  NS  equations  written  in 
generalized  coordinates,  equation  (25),  is  written  as  follows 


Q*.  ♦  E--  +  e;-  +  F-  ♦  F7-  =  0 


(  40  ) 


.  For  details  on  the  flux  splitting  used  in  this  algorithm 
see  reference  (7:1-9).  Equation  (40)  can  be  written  in  a 
factored  Beam-Warming  delta  form  as  follows 
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c5^"*’E  1  i  ,  j  -[  E  '  (  Qi+i/z,  j»Mi*i/a,  j  )-E'  (Qi-x/a,  j(Mx-i/a,  ;  )  ]  (42) 


Qi  +  x/a,  j-Qi  ,  j  +  cfoi  ,  j(  Qi.  j-Qi  -  i ,  j  )/2 


(  43  ) 


Qi  +  t/a,  j -Qi *i .  j  -  + 1 , j( Qi*a, j -Qx* i , j  )/2 


(  44  ) 


.  The  term  M  represents  geometric  terms  involved  in  the 
transformation  to  generalized  coordinates,  evaluated  at  the 
cell-interface  locations  where  the  flux  values  are  needed 
(7:9).  The  order  of  accuracy  of  the  approximation  is 
governed  by  the  value  of  the  switch  <p  :  second  order  accuracy 
for  cf)=  1,  first  order  for  (f)  =0  (7:4). 

The  terms  on  the  left  side  of  the  equality  in  equation 
(41)  are  solved  using  first  order  upwind  differencing,  which 
yields  a  block  tridiagonal  system  of  equations.  The  solution 
of  equation  (41)  is  achieved  by  factorization,  first  solving 
block-tr ldiagonal  equations  in  the  ^  direction  and  then  in 
the  Tj  direction  (7:9). 

The  fully  implicit  algorithm  was  written  assuming  no 
slip  conditions  at  the  body  surface,  adiabatic  wall,  zero 
pressure  gradient  normal  to  the  body  (19).  To  start  the 
solution  using  the  above  algorithm,  some  initial  conditions 
must  be  assumed.  Experience  gained  from  this  work  shows  a 
better  initialization  can  be  achieved,  i i  the  solution  for  a 
low  Mach  number,  is  used  to  initialize  the  flow  field  for  a 
high  Mach  number  solution.  If  the  final  solution  Mach  number 


15 


is  very  high,  repeated  cycles  from  lower  to  higher  Mach 
numbers  can  be  used.  The  details  of  this  initialization 
will  be  discussed  further  in  Chaper  IV. 

The  Upwind  Relaxation  Scheme 

When  solving  equation  <  30  )  with  subsonic  regions  in  the 
flow,  the  pressure  gradient  term  tangent  to  the  body  allows 
conditions  downstream  to  affect  conditions  upstream.  In 
subsonic  flow  regions  the  pressure  gradient  is  elliptic  in 
nature,  which  would  not  allow  the  solution  to  be  marched  in 
the  streamwise  direction.  To  overcome  this  difficulty, 
Vigneron,  Rakich  and  Tannehill  (27)  approximated  the 
streamwise  derivative  term  with  a  weighting  between  implicit 
and  explicit  differencing  that  depends  on  the  local  Mach 
number.  The  splitting  of  the  streamwise  flux  term  is  given 
as  follows 

E  =  E"  +  P  ( 45 


where 


E*  =  (  p  u,  pu^+COp,  p  uv,  (e*+p)u]T 
P  -  I  0,  (  1-  W>p,  0,  0  ]T 


(  46  ) 
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.  This  transformation  simplifies  the  (  E"  )  flux 

vector , equation  (31),  in  generalized  coordinates  to  the 

following 

E*  -  <  ^»EW+  £„P  1/J  (  58  ) 

,  and  equation  (30)  is  now  written  as  follows 

Q*  *  +  E*  +  F"  =  0  (59) 

Equation  (59)  is  a  system  of  nonlinear  differential 
equations  that  must  be  linearized  for  the  algorithm  being 
developed  to  make  use  of  well  known  techniques  for  solving 
systems  of  linear  differential  equations.  An  implicit 
differencing  scheme  will  be  used  that  repeatedly  sweeps  over 
the  computational  space  until  steady  state  is  reached.  This 
algorithm  is  similar  to  the  work  of  Thomas  et.  al.  (25), 
using  first  order  backward  spatial  differencing  on  the  ( E* ) 
flux  vector,  first  order  forward  spatial  differencing  on  the 
(P>  vector  and  second  order  implicit  central  spatial 
differencing  on  the  (  E,  Ev"  ,Fi  ,Fv."  )  flux  vectors.  Using  first 
order  forward  spatial  differencing  on  the  (P)  vector  allows 
conditions  upstream  to  propagate  downstream  in  regions  of 
subsonic  or  reverse  flow.  The  algorithm  finds  the  difference 
in  the  (Q)  vector  (Z^iQ)  such  that 


where  n  represents  the  time  level  of  the  solution.  Steady 
state  solutions  are  achieved  when  Aq  goes  to  zero. 

See  Appendix  B  for  details  on  the  linearization  and 
differencing  of  equation  (59),  which  yields  the  following 
equation 

(AlAQiw-*  +  IBlAQi.J  +  [ClAQi,J*i  - 
[  ^evj  iAq*-*.  j-  ^mE|*/j-t;wb?/j+7;wemV7/j- 
7)„Fi1/J+T7vF"v./?/J-  /J  (61) 

where  [A],[B]  and  [C]  are  matrices  defined  in  Appendix  B, 
equation  (b20),  subscripts  i  and  j  are  the  grid  node  indices 
for  the  ^  and  Tj  directions,  respectively. 

Viscous  Flux  Vector  Differencing 

The  viscous  flux  vectors  (F^"  )  contain  derivatives  with 
respect  to  7^,  and  are  being  differenced  with  repect  to  Tj  . 
Appendix  B  shows  the  differencing  technique  used  for  these 
terms,  and  the  following  is  sample  of  this  differencing 
technique.  The  second  element  in  the  (Fv"  )  vector  is 


differenced  as  follows 


[  (  4ya/3Re  )A,  4/Ll/3Re  )lt  j  ].5[  j+l, ** 

[Ui,i*i-Ui,j  )//\rTj  a- 

[  (  4ft /3Re  )i,  j+<  4ft/3Re)i.  j-A  ].5[  7}  a*/J  ]4f 

( ui.j-ui.j-i  i/At) m- 

[  (  2ft/3Re  )i  w^+<  2ft/3Re  )i,  j  1  .St  7) 

[Vi.JM-Vi.j  J/A^  a~ 

[  (  2/X/3Re  )i.  i+(  2ft/3Re  >i.  j-x  ].5[  T]H7)y/J  )i. 

(  vif  j-vi.  j-i  ]/At)  a  (62) 

.  It  should  be  noted  the  metrics  (  7j  Tj  )y,  and  the 

Jacobian  are  evaluated  at  j+1/2  and  j-1/2,  where  as  the 
coefficient  of  viscosity  and  Reynolds  number  are  determined 
by  taking  their  average  values  between  adjacent  grid  nodal 
points.  The  metrics  and  Jacobian  are  evaluated  at  half  nodal 
points . 

Three  grid  nodal  points  are  required  to  implicitly 
difference  the  viscous  terms  using  this  algorithm,  which 
results  in  solving  a  block  tridiagonal  matrix  to  get  the 
solution.  The  block  tridiagonal  solver  used  in  this 
algorithm  was  written  by  Carlson  (11). 

To  make  the  convergence  faster,  a  variable  time  step 
over  the  domain  is  implemented  similar  to  the  work  of  Halim 
and  Ghia  (13).  Such  techniques  are  always  used  to  enhance 
the  convergence  if  a  transient  solution  is  not  of  interest. 


IV  Results  and  Discussions 


Application  of  the  algorithms  described  in  Chapter  III 
to  a  series  of  high  speed  viscous  flow  problems  are  shown  in 
this  chapter.  The  cases  studied  using  the  fully  implicit 
scheme  are, 

-  The  blunt  body  flow 

Supersonic  viscous  flow  past  a  circular  wedge 
.  The  cases  studied  using  the  relaxation  scheme  are, 

Couette  flow  problem 
Supersonic  flow  over  flat  plate 
.  A  brief  description  of  each  problem  will  be  followed  by 
the  results  and  discussion  of  their  significance. 


The  Blunt  Body  Flow 

The  fully  implicit  algorithm  is  used  to  compute  the  two 
dimensional  flow  fields  surrounding  a  circular  cylinder. 
Although  the  present  results  are  two  dimensional,  qualitative 
comparisons  still  can  be  made  with  the  Edney  experiments  (9), 
especially  near  the  bow  shock  where  the  flow  is  locally 
two-dimensional . 

The  freestream  conditions  for  the  flow  are 

Re..  =  10000 

Pr  =0 . 72  /  =1.4 

p00=14 . 93  N/m* 

D-0. 3048m  T*=5S5.6  K 

The  grid  shown  in  Figure  1  was  used  for  this  computation 


M  =4.6 
oo 


T^_  =166.7  K 

OO 


and  has  81  grid  points  in  the  direction  tangent  to  the  body 


and  71  grid  points  in  the  direction  normal  to  the  body. 


The  entire  flow  field  was  initialized  at  the  freestream 
conditions.  It  was  difficult  to  get  the  solution  for  the 
blunt  body  problem  at  the  freestream  Mach  number  ( M^-4.6  ). 

To  overcome  this  difficulty  the  solution  was  obtained  for  a 
lower  Mach  number  and  stepped-up  to  its  final  value  over  1100 
iterations.  The  boundary  conditions  included  no  slip 
velocity  at  the  body  surface  and  no  pressure  gradient  normal 
to  the  body.  The  code  was  modified  slightly  to  be  able  to 
solve  for  constant  wall  temperature. 

Figure  2  shows  the  distribution  of  the 
nondimensional  wall  pressure  < Pw/P*%»* >  compared  to  the 
experimental  data  of  reference  (9).  Figure  3  shows  the 
distribution  of  the  nondimensional  heat  flux  (QM/Qata9  )  along 
with  the  boundary  layer  solution  of  reference  (26).  Plots  of 
Mach  contours,  velocity  vectors,  density  contours  and 


pressure  contours  are  shown  in  Figures  4-7. 
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Figure  2.  Comparison  of  Wall  Pressure 
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Figure  3.  Comparison  of  Heat  Flux 


The  comparison  of  the  wall  pressure  in  figure  2  and  the 
heat  transfer  in  figure  3  show  good  agreement  between  the 
fully  implicit  algorithm  and  the  available  experimental  and 
numerical  data.  That  gives  confidence  in  the  algorithm  to 
study  the  supersonic  flow  past  the  blunt  leading  edge 
geometry  of  reference  (10). 


Supersonic  Flow  Past  Circular  Wedge 


The  grid  used  for  the  initial  computation  has  81  grid 
points  in  the  direction  tangent  to  the  body  and  71  grid 
points  in  the  direction  normal  to  the  body.  The  solution  for 
the  initial  computation  was  then  used  as  the  initial 
condition  for  the  flow  over  a  modified  grid.  The  modified 
grid  clustered  points  near  the  body  surface  to  improve  the 
resolution  of  the  viscous  layer.  The  modified  grid  shown  in 
Figure  8  has  81  grid  points  tangent  to  the  body  and  62  grid 
points  normal  to  the  body.  The  grid  points  normal  to  the 
body  have  the  following  distribution;  32  evenly  spaced  grid 


points  in  the  viscous  layer  where  only  8  grid  points  had  been 


used  in  the  unmodified  grid  and  the  remaining  30  grid  points 
were  evenly  spaced  normal  to  the  body  over  the  rest  of  the 
grid. 

The  entire  flow  field  was  intiallized  at  the  freestream 
conditions.  It  was  difficult  to  get  the  solution  for  the 
circular  wedge  problem  at  the  freestream  Mach  number 
^^=6.57).  To  overcome  this  difficulty  the  solution  was 
obtained  for  a  lower  Mach  number  and  stepped-up  to  its  final 
value  over  1100  iterations.  No  slip  velocity  are  forced  at 
the  surface  of  the  body,  along  with  adiabatic  wall  and  no 
pressure  gradient  normal  to  the  body  surface.  Symmetry 
conditions  were  imposed  along  the  line  of  symmetry. 

Since  the  available  results  for  this  geometry  are  only 
from  an  inviscid  solution,  detailed  comparison  will  not  be 
possible  however  qualitative  preditions  can  be  seen.  For 
example  figure  9  shows  the  density  contours  using  the  present 
solution  compared  to  those  of  reference  (10). 

Mach  contours  are  shown  in  figure  10,  the  wall  pressure 
is  shown  in  figure  11,  and  the  velocity  profiles  at  the  grid 
location  X/R=l  are  shown  in  figures  12  and  13. 


The  next  set  of  results  are  obtained  using  the 
relaxation  algorithm  developed  in  the  present  study. 

One  Dimens ioal  Couette  Flow 

Couette  flow  is  known  as  the  flow  between  two  infinite 
parallel  plates  separated  by  a  distance  (H>.  To  study  this 
flow  the  relaxation  algorithm  was  solved  for  four  pressure 
gradients  in  the  direction  tangent  to  the  plates.  This  test 
case  was  run  to  check  the  formulation  and  the  numerical 
procedures  of  the  algorithm. 

The  xower  plate  was  held  fixed  while  the  upper  plate 
moved  at  Mach  0.09.  The  Reynolds  number  for  this  flow  is 
6.2,  and  the  four  pressure  gradients  analysed  such  that  B- 
-2, 2, 4, 7;  where  B"-<  Re«  )<  p  )*/[  4(  M^  )  )  (28:113).  To  convert  B 
to  correspond  to  nondimensionalization  used  in  this 
algorithm,  B*  is  defined  where  B’=4(M)(B). 

The  computational  grid  for  this  test  case  used  12  grid 
points  in  the  direction  normal  to  the  plates  with  y/H=l,  and 
the  grid  points  in  the  direction  tangent  to  the  plates  were 
varied  between  10,20  and  30  with  x/H=l. 

Initially  the  entire  flow  was  given  the  same  velocity 
profile  and  the  desired  pressure  gradient.  No  slip  velocity 
conditions  are  forced  at  the  upper  and  lower  plates,  while 
the  Initial  conditions  at  grid  location  i=l  were  not  allowed 
to  change  during  the  computation  (i  is  the  grid  position 
index  in  the  direction  tangent  to  the  plate). 

After  computing  the  solution  for  the  four  pressure 


gradients  for  each  of  the  three  different  x  grid  spacings, 
the  solutions  at  grid  point,  y/H-.8  were  compared  (see  figure 
14).  As  the  spacing  between  grid  points  decreases,  the 
solution  calculated  by  the  algorithm  should  more  closely 
approach  the  exact  solution.  Figure  14,  shows  that  for  the 
three  x  grid  spacings  chosen  the  solution  agreed  exactly  with 
the  analytic  solution  (DELTA  X  =  0.0).  As  a  side  note, 
although  the  NS  equations  in  general  do  not  have  an  exact 
solution,  the  exact  solution  at  specific  points  can  be 
obtained  numerically  using  the  accuracy  study  as  shown  in 
figure  14. 

A  comparison  between  the  velocity  profiles  calculated  by 
the  algorithm  and  the  profiles  calculated  using  the 
analytical  solution  are  presented  in  figure  15. 
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Figure  14.  Accuracy  study  of  Couette  flow  for 
different  pressure  gradients  (B>. 


The  second  test  problem  for  the  relaxation  algorithm  is 


solving  the  supersonic  flow  over  a  flat  plate. 

The  freestream  conditions  for  the  flow  are 

=  2 

OO 

Re  /L**  =  1.65X10*/m 
oo 

L-  =  lm 

(Too  >"  =  <T->-'  *  221.6  K 

Pr  =  .72 
Y  *  1.4 

The  grid  used  41  grid  points  normal  to  the  plate 
surface,  with  a  DELTA( y  )= . 1524X10-3m,  which  produced  a 
constant  grid  height  of  .61X10-am.  The  grid  used  52  grid 
points  tangent  to  the  plate  surface,  with  the  initial  grid 
point  i  =  l  relating  to  the  x  distance  of  0.305m  from  the 
leading  edge  of  the  plate.  The  1=52  grid  point  relating  to 
the  x  distance  of  0.915m  from  the  leading  edge  of  the  plate, 
requiring  DELTA< x  )=0 . 0 122m . 

The  entire  flow  field  was  mitiallized  with  the  velocity 
and  temperature  profiles  calculated  at  plate  position  0.305m 
by  a  boundary  layer  code  written  by  Cebeci  and  Bradshaw  (12). 
No  slip  velocity  conditions  were  forced  at  the  plate  surface, 
and  the  normal  velocity  component  was  allowed  to  propagate 
out  the  top  of  the  grid.  The  initial  flow  conditions  at  grid 
location  i=l  were  not  allowed  to  change  during  the 
computat ions . 

The  tangential  velocity  and  temperature  profiles  at 


plate  position  0.915m  are  compared  to  the  Cebeci  boundary 
layer  solution  at  the  same  location  (see  figure  16  and  figure 
17  ). 

As  seen,  relaxation  algorithms  would  allow  someone  to 
take  larger  step  sizes  in  the  streamwise  direction  compared 
to  space  marching  procedures.  For  example,  the  same  problem 
is  studied  by  Huband  (14)  using  a  space  marching  procedure. 
The  number  of  points  used  in  the  streamwise  direction  are  612 
compared  to  51  used  in  the  present  study,  however  the  CPU 
time  taken  to  get  the  final  solution  using  the  present  method 
is  much  larger  compared  to  the  space  marching  procedure. 
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V  Conclusions  and  Recommendations 

The  solutions  obtained  using  the  implicit  Navier-Stokes 
algorithm,  compared  favorably  with  experimental  data.  It  is 
also  recommended  that  this  algorithm  be  modified  to  solve 
real  gas  problems. 

The  Approximate  Navier-Stokes  Algorithm  developed  in 
this  study  also  agreed  well  with  analytic  solutions  and  other 


numeric  solutions.  The  time  required  to  achieve  these 
solutions  could  be  reduced  by  updating  the  Jacobian  metrices 


To  transform  the  governing  equations  into  generalized 


coordinates  the  chain  rule  is  used  to  determine  the 
derivatives  of  the  general  coordiante  variables  with  respect 
to  the  variables  in  physical  space.  Let  the  following 
general  transforms  be  used 

t  *t  ( al  ) 

^  x,y  )  (  a2 ) 

'V  -  77<xty>  (  a3  ) 


with 


(  )te  =  <  >*  (  a4  ) 

(  >*  =  (  ^  )„<  +(  Tj  )„<  (  a5  ) 

(  )y  =  <  ^  M  +  <  'T?  M  >rj  (  a6  ) 

.  The  unit  area  in  physical  space  is  related  to  the  unit 
area  in  the  general  coordinate  space  by  the  Jacobian  of 
transformation  (J).  J  is  the  determinate  of  the  matrix 
formed  by  the  metrics  (  ^  >*,(  Tj  )„,(  ^  )y  and  (  '7]  )y ,  which  is 


r’nx  =  (  <  2  +277  )-(  ^  yv§ +07yV?  >  12(a/3Re  (  a23  ) 

T'*y  *  (  ^>u^+77»,u^+  ^xV^+7 )/il/Re  (  a24  ) 

r'yy  -  [  (2  ^yv^  +2  7]yV^  )-(  f  *u^  +  77*u,j  >  \2{J,/3Re  (  a25  ) 

q ' m  -  «AU  +  7)>.T?  ]/[  RePr (  /  -1  >  ]  (  a26  ) 

q'y  =  -yLU  ^yT^7;yT?  ]/[  RePrC  7 -1  )  ]  (  a27  ) 

where 

Re  *  (a28) 

Pr  =  cP*/4*/ k4*  (  a29  ) 


cP  -  /R/<  /  -1  ) 


(  a 30  ) 


Appendix  B:  Details  about  the  Upwind  Relaxation  Scheme 
Linearization  of  the  ANS  Equations 

To  linearize  equation  (59)  of  the  text 

Q%  +  B~  +  F^  *  0  ( bl  ) 

use  the  following  Taylor  expansions 


E*"*1  =  E«r.  +  A^"B«* 

=  +  A^"Q$E-« 

-  E**"  +  AQ"E“^  (  b2  ) 

where  ( Q  can  be  written  as  Aq/A f.  and  n  denotes  the  time 
level  of  the  calculation.  In  like  manner  the  linearization 
of  the  other  terms  are  given  as 


E"* 1  “  E"  +  AQ”E«  (  b3  ) 

Ev"-1  =  +  AQnEv«  (  b4  ) 

Ft"*1  =  Fi"  +AQ"Fia,  (b5) 

Fv"*‘  =  Fv"  +  AQnFv«  (  b6  ) 


.  Substituting  equations  (b2,b3,b4,b5  and  b6  )  into  equations 
(58  and  32)  which  are  then  substituted  into  equation  (  bl  ) 


yeilding 


Ao/At/j  +  t  Aq  f-E%/j  ij  +  (Aq^kE^/j  in  - 

|Aq/7]-Emv«/J  1,  +(AQ^vFi«/J  1^-[A<3'??xF”v*/J  1,  - 
-  ^wE»/J-/T7ME1/J+7}«E"^^/J-07yFi?/J  +  07>F”v^/J-  fxP^/J  <  b7  ) 

where  only  terms  containing  A©  are  retained  on  the  left  hand 
side  ( LHS  )  of  the  equality,  and  these  terms  will  be 
differenced  implicitly  at  the  < n+1  )  time  level  to  produce  a 
block  tridiagonal  matrix  multiplied  by  the  A®  vector  at  the 
(n)  time  level.  The  terms  on  the  right  hand  side  ( RHS >  of 
the  equality  are  all  known  at  the  (n)  time  level. 

Finite  Difference  Equations 

The  ANS  equations,  written  in  vector  form  in  equation 
<b7),  are  differenced  term  by  term  as  follows 

AQ/At/J  =  (  [  I  J/At/J  )i,  .  J  <b0) 

|Aq^hE%/J  !?  =  t(§xE*»,/J)i.,JA(2i.J- 

<  £;„e%/j  )*_» ,  jAQi-» .  j  i/Af  (  by  > 

i  Aq^kEq/j  )?  =  [  (  nrj.E^/J  )i,  j**A<3i.  j-i- 

(  7]KEm/J  >i,  J-Ai  •  J-i  l/2AT 


(  blO  ) 


F”v  =  [0,  r’-y,  rVy,  ur'xy+vr'^r  <  b30  > 

r'yy  =  4  fJLV)v<,  Qs/Qj.  )fj/3Re-2  Qa/Q»  )rj/3Re  (b31) 

qx  =  fA-Y T/^l  <Q4/Qi  )-(  Q**+Q3*  >/2Q»  l^/RePr  (  b32  ) 

with  these  definitions  the  Jacobian  matrices  can  be  formed  as 
shown  on  the  following  pages. 
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s,V\^,‘n'  * 


11 

12 

13 

14 

21 

22 

23 

24 

31 

32 

33 

34 

41 

42 

43 

44 

11  =  0 
12  =  1 
13-0 
14  =  0 

21  «  -Qaa/Qxa  ♦  (  / -1  )<Q»a+Q3a  )/2Qi* 

22  -  2Q»/Q i  -  (/-l)Qa/Qi 

23  -  -(  Y -1  )Qa/Qi 

24  -  (  /-I  ) 

31  -  -QaQa/Qxa 

32  *  Q»/Qx 

33  -  Qa/Q. 

34-0 

41  =  -Qa/Qx*{Q4><  X-l  >1  Q4-<  Qa*+Qsa  >/2Q» 

+  Qa/Qi  [  (  X-l  X  Qaa+Qsa  >2Qi2  ] 

42  =  1/Qi{  Q*+<  X-l  >(  Q4-<  Qaa+Qsa  >/2Q»  ]  > 

-  </-!)(  Qa/Qi  >a 


(  b34  ) 


43  -  -<  X -1 >QsQ*/Q. 

44  -  C  X >Qa/Q» 


<FV"  >*=  11 

12 

13 

14  ( b37  ) 

21 

22 

23 

24 

31 

32 

33 

34 

41 

42 

43 

44 

11  -  0 
12  =  0 
13*0 
14  *  0 

21  -  -(  fJL  X  7)  >„[  Qa/Q*a  ]^/<  Re  ) 

-  (  /X  X  7 )  M  Qa/Q  i a  ]#j/<  Re  ) 

22  =  (  /l  )(  7]  )yf  1/Q*  ],)/<  Re  ) 

23  *  (  /X  X  ^  )*[  1/Q*  ],}/<  Re  ) 

24-0 

31  -  -4(  fJ  x7))y[Qa/Qx*],j/3(Re) 

+  2(  /i  X  T;  )*[  Qa/Q  1  *  ]rj/3(  Re  ) 

32  =  -2<  /IX  7}  )x[  1/Q 1  }?/3(  Re  ) 

33  =  4(/iX7]  )y[  1/Q*  ]?/3(  Re  ) 

34*0 

41  -  -(/Ux7)MQa*/Q,*|(|/(Re> 

-4(  /IX  T;  )x[  Q,*/Q*»  J^/3(  Re  ) 

-(  fJL  X  )K{  (  Qa/Q  1  a  )[  Q3/Q 1  1/|  +  <  Qa/Q»  M  Qa/Q*  ]r|  ) 

/(  Re  ) 

+  <  H  X  39  )*{  (Qa/Q*2  >f  Qa/Q*  ir,+(Q3/Q»  )[  Qa/Q*  a  ] ^  > 

/3(  Re  ) 

t(  U  )(  y  )(  71  :„[  -<  Q«/Q»  *  >♦<  Qaa+Qaa  >/Q*3  )«/<  Re  )(  Pr  > 


(  fJ.  X  7]  )y[  Qa/Qi a  }rj  /<  Re  ) 

+(^)(7)  >*[  Q3/Q1  ]?/(  Re  )Qi 
-2(  /i)(7)  )xQ3[  1/Q»  ]  ?/3(  Re  )Q» 
-(/iX/x'T)  >y[  Qa/Qi3  ]?/<  Re  )(  Pr  ) 
4<  yU.  )(7)  )y[  Qs/Qi  a  ]»j/3(  Re  ) 

+(  yU.  X  7}  )«Qa[  1/Qi  ]^/<  Re  >Qi 
-2(  yU  X  7?  >M[  Qa/Qi  ]q/3<  Re  )Q» 

-<  yU  )<  /  )(  7)  >„[  Qs/Qi3  ]/j/(  Re  )(  Pr  ) 
(  fJ.  X  /  X  'Tj  )y[  1/Qi  )^/<  Re  X  Pr  ) 
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